module GOLD_PointAccel
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of GOLD.                                        *
!*                                                                     *
!* GOLD is free software; you can redistribute it and/or modify it and *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* GOLD is distributed in the hope that it will be useful, but WITHOUT *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

!***********************************************************************
!*                                                                     *
!*     The two subroutines in this file write out all of the terms     *
!*  in the u- or v-momentum balance at a given point.  Usually         *
!*  these subroutines are called after the velocities exceed some      *
!*  threshold, in order to determine which term is culpable.           *
!*  often this is done for debugging purposes.                         *
!*                                                                     *
!*  Macros written all in capital letters are defined in GOLD_memory.h *
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q, f                                     *
!*    j+1  > o > o >   At ^:  v, PFv, CAv, vh, diffv, vbt, vhtr        *
!*    j    x ^ x ^ x   At >:  u, PFu, CAu, uh, diffu, ubt, uhtr        *
!*    j    > o > o >   At o:  h, D, s, tr, T, S                        *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

use GOLD_diag_mediator, only : diag_ptrs
use GOLD_domains, only : pe_here
use GOLD_error_handler, only : GOLD_error, NOTE
use GOLD_file_parser, only : get_param, log_param, log_version, param_file_type
use GOLD_grid, only : ocean_grid_type
use GOLD_io, only : create_file, open_file, reopen_file
use GOLD_io, only : APPEND_FILE, ASCII_FILE, MULTIPLE, SINGLE_FILE
use GOLD_time_manager, only : time_type, get_time, get_date, set_date, operator(-)
use GOLD_variables, only : directories, ocean_internal_state

implicit none ; private

#include <GOLD_memory.h>

public write_u_accel, write_v_accel, PointAccel_init

type, public :: PointAccel_CS ; private
  logical :: split          ! If true, use the split time stepping scheme.
  character(len = 200) :: u_trunc_file ! The complete path to files in which a
  character(len = 200) :: v_trunc_file ! column's worth of accelerations are
                                       ! written if velocity truncations occur.
  integer :: u_file, v_file ! The unit numbers for opened u- or v- truncation
                            ! files, or -1 if they have not yet been opened.
  integer :: cols_written   ! The number of columns whose output has been
                            ! written by this PE during the current run.
  integer :: max_writes     ! The maximum number of times any PE can write out
                            ! a column's worth of accelerations during a run.
  type(time_type), pointer :: Time ! A pointer to the ocean model's clock.
  type(diag_ptrs), pointer :: diag ! A pointer to a structure of shareable
                            ! ocean diagnostic fields.
! The following are pointers to many of the state variables and accelerations
! that are used to step the physical model forward.  They all use the same
! names as the variables they point to in GOLD.F90
  real, pointer, dimension(:,:,:,:) :: &
    u => NULL(), v => NULL()     ! Velocities in m s-1.
  real, pointer, dimension(:,:,:) :: &
    u_av => NULL(), v_av => NULL(), & ! Time average velocities in m s-1.
    u_prev => NULL(), v_prev => NULL(), & ! Previous velocities in m s-1.
    uh => NULL(), vh => NULL(), &   ! Volume fluxes in m3 s-1.
    CAu => NULL(), CAv => NULL(), & ! Coriolis accelerations in m s-2.
    PFu => NULL(), PFv => NULL(), & ! Pressure force accelerations in m s-2.
    diffu => NULL(), diffv => NULL(), & ! Horizontal viscous accelerations.
    T => NULL(), S => NULL(), &     ! Temperature and salinity in C and psu.
    pbce => NULL(), &               ! pbce times eta gives the baroclinic
                                    ! pressure anomaly in each layer due to
                                    ! free surface height anomalies.
                                    ! pbce has units of m s-2.
    u_accel_bt => NULL(), &         ! Barotropic acclerations in m s-2.
    v_accel_bt => NULL()

end type PointAccel_CS

contains

subroutine write_u_accel(I, j, um, hin, dt, G, CS, maxvel, minvel, str, a, hv)
  integer,                                intent(in) :: I, j
  real, dimension(NXMEMQ_,NYMEM_,NZ_),    intent(in) :: um
  real, dimension(NXMEM_,NYMEM_,NZ_),     intent(in) :: hin
  real,                                   intent(in) :: dt
  type(ocean_grid_type),                  intent(in) :: G
  type(PointAccel_CS),                    pointer    :: CS
  real,                                   intent(in) :: maxvel, minvel
  real, optional,                         intent(in) :: str
  real, dimension(NXMEMQ_,NZ_), optional, intent(in) :: a, hv
! This subroutine writes to an output file all of the accelerations
! that have been applied to a column of zonal velocities over the
! previous timestep.  This subroutine is called from vertvisc.

! Arguments: I - The zonal index of the column to be documented.
!  (in)      j - The meridional index of the column to be documented.
!  (in)      um - The new zonal velocity, in m s-1.
!  (in)      hin - The layer thickness, in m.
!  (in)      dt - The model's dynamics time step.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 PointAccel_init.
!  (in)      str - The surface wind stress integrated over a time
!                  step, in m2 s-1.
!  (in)      a - The layer coupling coefficients from vertvisc, m.
!  (in)      hv - The layer thicknesses at velocity grid points, from
!                 vertvisc, in m.

  real, pointer, dimension(:,:,:) :: &
    u_prev, &     ! The previous zonal velocity, in m s-1.
    v_prev        ! The previous meridional velocity, in m s-1.
  real    :: f_eff, CFL
  real    :: truncvel, du
  real    :: Inorm(SZK_(G))
  real    :: e(SZK_(G)+1)
  integer :: yr, mo, day, hr, minute, sec, yearday
  integer :: k, m, mp, ks, ke
  integer :: nz
  logical :: do_k(SZK_(G)+1)
  integer :: file

!  if (.not.associated(CS)) return
  nz = G%ke
  if (CS%cols_written < CS%max_writes) then
    CS%cols_written = CS%cols_written + 1

    ks = 1 ; ke = nz ; m = 1
    do_k(:) = .false.

  ! Open up the file for output if this is the first call.
    if (CS%u_file < 0) then
      if (len_trim(CS%u_trunc_file) < 1) return
      call open_file(CS%u_file, trim(CS%u_trunc_file), action=APPEND_FILE, &
                     form=ASCII_FILE, threading=MULTIPLE, fileset=SINGLE_FILE)
      if (CS%u_file < 0) then
        call GOLD_error(NOTE, 'Unable to open file '//trim(CS%u_trunc_file)//'.')
        return
      endif
    endif
    file = CS%u_file

    if ( .not.(associated(CS%u_prev) .and. associated(CS%v_prev)) ) then
      ! Determine the index in u of um.
      do k=1,nz
        if (CS%u(I,j,k,1) /= CS%u(I,j,k,2)) then
          if (um(I,j,k) == CS%u(I,j,k,1)) then ; m = 1
          elseif (um(I,j,k) == CS%u(I,j,k,2)) then ; m = 2
          else
            write(file,'("Mismatched u velocities in GOLD_pointaccel")') ; exit
          endif
          exit
        endif
      enddo
      mp = 3 - m
    endif

    if (associated(CS%u_prev)) then
      u_prev => CS%u_prev
    else
      u_prev => CS%u(:,:,:,mp)
    endif
    if (associated(CS%v_prev)) then
      v_prev => CS%v_prev
    else
      v_prev => CS%v(:,:,:,mp)
    endif

  ! Determine which layers to write out accelerations for.
    do k=1,nz
      if (((max(CS%u_av(I,j,k),um(I,j,k)) >= maxvel) .or. &
           (min(CS%u_av(I,j,k),um(I,j,k)) <= minvel)) .and. &
          ((hin(i,j,k) + hin(i+1,j,k)) > 3.0*G%Angstrom)) exit
    enddo
    ks = k
    do k=nz,1,-1
      if (((max(CS%u_av(I,j,k), um(I,j,k)) >= maxvel) .or. &
           (min(CS%u_av(I,j,k), um(I,j,k)) <= minvel)) .and. &
          ((hin(i,j,k) + hin(i+1,j,k)) > 3.0*G%Angstrom)) exit
    enddo
    ke = k
    if (ke < ks) then
      ks = 1; ke = nz; write(file,'("U: Unable to set ks & ke.")')
    endif

    call get_date(CS%Time, yr, mo, day, hr, minute, sec)
    call get_time((CS%Time - set_date(yr, 1, 1, 0, 0, 0)), sec, yearday)
    write (file,'(/,"--------------------------")')
    write (file,'(/,"Time ",i5,i4,F6.2," U-velocity violation at ",I4,": ",2(I3), &
        & " (",F7.2," E "F7.2," N) Layers ",I3," to ",I3,". dt = ",1PG10.4)') &
        yr, yearday, (REAL(sec)/3600.0), pe_here(), I, j, &
        G%geolonu(I,j), G%geolatu(I,j), ks, ke, dt

    if (ks <= G%nk_rho_varies) ks = 1
    do k=ks,ke
      if ((hin(i,j,k) + hin(i+1,j,k)) > 3.0*G%Angstrom) do_k(k) = .true.
    enddo

    write(file,'(/,"Layers:",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(I10," ",$)') (k); enddo
    write(file,'(/,"u(m):  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (um(I,j,k)); enddo
    write(file,'(/,"u(mp): ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (u_prev(I,j,k)); enddo
    write(file,'(/,"u(3):  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (CS%u_av(I,j,k)); enddo

    write(file,'(/,"CFL u: ",$)')
    do k=ks,ke ; if (do_k(k)) then
      CFL = abs(um(I,j,k)) * dt * G%dy_u(I,j)
      if (um(I,j,k) < 0.0) then ; CFL = CFL * G%IDXDYh(i+1,j)
      else ; CFL = CFL * G%IDXDYh(i,j) ; endif
      write(file,'(ES10.3," ",$)') CFL
    endif ; enddo
    write(file,'(/,"CFL0 u:",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    abs(um(I,j,k)) * dt * G%IDXu(I,j) ; enddo

    write(file,'(/,"du:    ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    ((um(I,j,k)-u_prev(I,j,k))); enddo
    write(file,'(/,"CAu:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (dt*CS%CAu(I,j,k)); enddo
    write(file,'(/,"PFu:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (dt*CS%PFu(I,j,k)); enddo
    if (ASSOCIATED(CS%diag%PFu_bc)) then
      write(file,'(/,"PFubc: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (dt*CS%diag%PFu_bc(I,j,k)); enddo
    endif
    write(file,'(/,"diffu: ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (dt*CS%diffu(I,j,k)); enddo
    if (ASSOCIATED(CS%diag%Ah_h)) then
      write(file,'(/,"Ah_h+:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%Ah_h(i,j,k)); enddo
      write(file,'(/,"Ah_h-:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%Ah_h(i-1,j,k)); enddo
    endif
    if (ASSOCIATED(CS%diag%Ah_q)) then
      write(file,'(/,"Ah_q+:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%Ah_q(I,J+1,k)); enddo
      write(file,'(/,"Ah_q-:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%Ah_q(I,J,k)); enddo
    endif

    if (ASSOCIATED(CS%diag%Kh_h)) then
      write(file,'(/,"Kh_h+:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%Kh_h(i,j,k)); enddo
      write(file,'(/,"Kh_h-:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%Kh_h(i-1,j,k)); enddo
    endif
    if (ASSOCIATED(CS%diag%Kh_q)) then
      write(file,'(/,"Kh_q+:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%Kh_q(I,J+1,k)); enddo

      write(file,'(/,"Kh_q-:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%Kh_q(I,J,k)); enddo
    endif

    if (ASSOCIATED(CS%diag%gradKEu)) then
      write(file,'(/,"KEu:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (dt*CS%diag%gradKEu(I,j,k)); enddo
    endif
    if (ASSOCIATED(CS%diag%rv_x_v)) then
      write(file,'(/,"Coru:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
          dt*(CS%CAu(I,j,k)-CS%diag%rv_x_v(I,j,k)); enddo
    endif
    if (ASSOCIATED(CS%diag%du_dt_visc)) then
      write(file,'(/,"ubv:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
          (um(I,j,k)-dt*CS%diag%du_dt_visc(I,j,k)); enddo
      write(file,'(/,"duv:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (dt*CS%diag%du_dt_visc(I,j,k)); enddo
    endif
    if (ASSOCIATED(CS%diag%du_adj)) then
      write(file,'(/,"duadj: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%du_adj(I,j,k)); enddo
    endif
    if (ASSOCIATED(CS%diag%du_adj2)) then
      write(file,'(/,"duadj2: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%du_adj2(I,j,k)); enddo
    endif
    if (present(a)) then
      write(file,'(/,"a:     ",$)')
      do k=ks,ke+1 ; if (do_k(k)) write(file,'(ES10.3," ",$)') a(I,k); enddo
    endif
    if (present(hv)) then
      write(file,'(/,"hvel:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') hv(I,k); enddo
    endif
    write(file,'(/,"Stress:  ",ES10.3)') str

    if (CS%split) then
      write(file,'("dubt:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (dt*CS%u_accel_bt(I,j,k)) ; enddo
      write(file,'(/)')
      if (ASSOCIATED(CS%diag%PFu_bt)) &
        write(file,'("PFubt: ",ES10.3)') (dt*CS%diag%PFu_bt(I,j))
      if (ASSOCIATED(CS%diag%Coru_bt)) &
        write(file,'("Corubt:",ES10.3)') (dt*CS%diag%Coru_bt(I,j))
      if (ASSOCIATED(CS%diag%Nonlnu_bt)) &
        write(file,'("NLubt: ",ES10.3)') (dt*CS%diag%Nonlnu_bt(I,j))
    endif

    write(file,'(/,"h--:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (hin(i,j-1,k)); enddo
    write(file,'(/,"h+-:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (hin(i+1,j-1,k)); enddo
    write(file,'(/,"h-0:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (hin(i,j,k)); enddo
    write(file,'(/,"h+0:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (hin(i+1,j,k)); enddo
    write(file,'(/,"h-+:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (hin(i,j+1,k)); enddo
    write(file,'(/,"h++:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (hin(i+1,j+1,k)); enddo


    e(nz+1) = -G%D(i,j)
    do k=nz,1,-1 ; e(K) = e(K+1) + hin(i,j,k) ; enddo
    write(file,'(/,"e-:    ",$)')
    write(file,'(ES10.3," ",$)') e(ks)
    do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ",$)') e(K); enddo

    e(nz+1) = -G%D(i+1,j)
    do k=nz,1,-1 ; e(K) = e(K+1) + hin(i+1,j,k) ; enddo
    write(file,'(/,"e+:    ",$)')
    write(file,'(ES10.3," ",$)') e(ks)
    do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ",$)') e(K) ; enddo
    if (ASSOCIATED(CS%T)) then
      write(file,'(/,"T-:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') CS%T(i,j,k); enddo
      write(file,'(/,"T+:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') CS%T(i+1,j,k); enddo
    endif
    if (ASSOCIATED(CS%S)) then
      write(file,'(/,"S-:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') CS%S(i,j,k); enddo
      write(file,'(/,"S+:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') CS%S(i+1,j,k); enddo
    endif

    write(file,'(/,"v--:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (v_prev(i,J-1,k)); enddo
    write(file,'(/,"v-+:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (v_prev(i,J,k)); enddo
    write(file,'(/,"v+-:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (v_prev(i+1,J-1,k)); enddo
    write(file,'(/,"v++:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (v_prev(i+1,J,k)); enddo


    write(file,'(/,"vh--:  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (CS%vh(i,J-1,k)*G%IDXv(i,J-1)); enddo
    if (ASSOCIATED(CS%diag%vh_max)) then
      write(file,'(/," vh^--: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%vh_max(i,J-1,k)*G%IDXv(i,J-1)); enddo
    endif
    if (ASSOCIATED(CS%diag%vh_min)) then
      write(file,'(/," vh_--: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%vh_min(i,J-1,k)*G%IDXv(i,J-1)); enddo
    endif
    if (ASSOCIATED(CS%diag%vh_lay)) then
      write(file,'(/," vhL--: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%vh_lay(i,J-1,k)*G%IDXv(i,J-1)); enddo
    endif
    if (ASSOCIATED(CS%diag%vh_cent)) then
      write(file,'(/," vhC--: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%vh_cent(i,J-1,k)*G%IDXv(i,J-1)); enddo
    endif
    write(file,'(/," vhCp--:",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                        (0.5*CS%v_av(i,j-1,k)*(hin(i,j-1,k) + hin(i,j,k))); enddo

    write(file,'(/,"vh-+:  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (CS%vh(i,J,k)*G%IDXv(i,J)); enddo
    if (ASSOCIATED(CS%diag%vh_max)) then
      write(file,'(/," vh^-+: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%vh_max(i,J,k)*G%IDXv(i,J)); enddo
    endif
    if (ASSOCIATED(CS%diag%vh_min)) then
      write(file,'(/," vh_-+: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%vh_min(i,J,k)*G%IDXv(i,J)); enddo
    endif
    if (ASSOCIATED(CS%diag%vh_lay)) then
      write(file,'(/," vhL-+: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%vh_lay(i,J,k)*G%IDXv(i,J)); enddo
    endif
    if (ASSOCIATED(CS%diag%vh_cent)) then
      write(file,'(/," vhC-+: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%vh_cent(i,J,k)*G%IDXv(i,J)); enddo
    endif
    write(file,'(/," vhCp-+:",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                        (0.5*CS%v_av(i,J,k)*(hin(i,j,k) + hin(i,j+1,k))); enddo

    write(file,'(/,"vh+-:  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%vh(i+1,J-1,k)*G%IDXv(i+1,J-1)); enddo
    if (ASSOCIATED(CS%diag%vh_max)) then
      write(file,'(/," vh^+-: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                          (CS%diag%vh_max(i+1,J-1,k)*G%IDXv(i+1,J-1)); enddo
    endif
    if (ASSOCIATED(CS%diag%vh_min)) then
      write(file,'(/," vh_+-: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                          (CS%diag%vh_min(i+1,J-1,k)*G%IDXv(i+1,J-1)); enddo
    endif
    if (ASSOCIATED(CS%diag%vh_lay)) then
      write(file,'(/," vhL+-: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                          (CS%diag%vh_lay(i+1,J-1,k)*G%IDXv(i+1,J-1)); enddo
    endif
    if (ASSOCIATED(CS%diag%vh_cent)) then
      write(file,'(/," vhC+-: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                          (CS%diag%vh_cent(i+1,J-1,k)*G%IDXv(i+1,J-1)); enddo
    endif
    write(file,'(/," vhCp+-:",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                    (0.5*CS%v_av(i+1,J-1,k)*(hin(i+1,j-1,k) + hin(i+1,j,k))); enddo

    write(file,'(/,"vh++:  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                          (CS%vh(i+1,J,k)*G%IDXv(i+1,J)); enddo
    if (ASSOCIATED(CS%diag%vh_max)) then
      write(file,'(/," vh^++: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                          (CS%diag%vh_max(i+1,J,k)*G%IDXv(i+1,J)); enddo
    endif
    if (ASSOCIATED(CS%diag%vh_min)) then
      write(file,'(/," vh_++: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                          (CS%diag%vh_min(i+1,J,k)*G%IDXv(i+1,J)); enddo
    endif
    if (ASSOCIATED(CS%diag%vh_lay)) then
      write(file,'(/," vhL++: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                          (CS%diag%vh_lay(i+1,J,k)*G%IDXv(i+1,J)); enddo
    endif
    if (ASSOCIATED(CS%diag%vh_cent)) then
      write(file,'(/," vhC++: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                          (CS%diag%vh_cent(i+1,J,k)*G%IDXv(i+1,J)); enddo
    endif
    write(file,'(/," vhCp++:",$)')
         do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                     (0.5*CS%v_av(i+1,J,k)*(hin(i+1,j,k) + hin(i+1,j+1,k))); enddo


    if (ASSOCIATED(CS%diag%q)) then
      write(file,'(/,"v-eff: ",$)')
      f_eff = G%f(I,J-1)
      if (abs(f_eff) < 1e-10) f_eff = 1e-10
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') 0.5*(CS%diag%q(I,J-1,k) * &
        (CS%vh(i,J-1,k) + CS%vh(i+1,J-1,k)) * G%IDXu(I,j) / f_eff); enddo

      write(file,'(/,"v+eff: ",$)')
      f_eff = G%f(I,J)
      if (abs(f_eff) < 1e-10) f_eff = 1e-10
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') 0.5*(CS%diag%q(I,J,k) * &
        (CS%vh(i+1,J,k) + CS%vh(i,J,k)) * G%IDXu(I,j) / f_eff); enddo
    endif

    write(file,'(/,"D:     ",2(ES10.3))') G%D(i,j),G%D(i+1,j)

  !  From here on, the normalized accelerations are written.

    do k=ks,ke
      du = um(I,j,k)-u_prev(I,j,k)
      if (abs(du) < 1.0e-6) du = 1.0e-6
      Inorm(k) = 1.0 / du
    enddo

    write(file,'(2/,"Norm:  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') (1.0/Inorm(k)); enddo

    write(file,'(/,"du:    ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                      ((um(I,j,k)-u_prev(I,j,k))*Inorm(k)); enddo

    write(file,'(/,"CAu:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                    (dt*CS%CAu(I,j,k)*Inorm(k)); enddo

    write(file,'(/,"PFu:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                    (dt*CS%PFu(I,j,k)*Inorm(k)); enddo

    if (ASSOCIATED(CS%diag%PFu_bc)) then
      write(file,'(/,"PFubc: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                      (dt*CS%diag%PFu_bc(I,j,k)*Inorm(k)); enddo
    endif

    write(file,'(/,"diffu: ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                    (dt*CS%diffu(I,j,k)*Inorm(k)); enddo

    if (ASSOCIATED(CS%diag%gradKEu)) then
      write(file,'(/,"KEu:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                      (dt*CS%diag%gradKEu(I,j,k)*Inorm(k)); enddo
    endif
    if (ASSOCIATED(CS%diag%rv_x_v)) then
      write(file,'(/,"Coru:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
          dt*(CS%CAu(I,j,k)-CS%diag%rv_x_v(I,j,k))*Inorm(k); enddo
    endif
    if (ASSOCIATED(CS%diag%du_dt_visc)) then
      write(file,'(/,"duv:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
          (dt*CS%diag%du_dt_visc(I,j,k))*Inorm(k); enddo
    endif
    if (ASSOCIATED(CS%diag%du_adj)) then
      write(file,'(/,"duadj: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
          (CS%diag%du_adj(I,j,k))*Inorm(k); enddo
    endif
    if (ASSOCIATED(CS%diag%du_adj2)) then
      write(file,'(/,"duadj2:",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
          (CS%diag%du_adj2(I,j,k))*Inorm(k); enddo
    endif
    if (CS%split) then
      write(file,'(/,"dubt:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                      (dt*CS%u_accel_bt(I,j,k)*Inorm(k)) ; enddo
    endif

    call flush(file)
  endif

end subroutine write_u_accel


subroutine write_v_accel(i, J, vm, hin, dt, G, CS, maxvel, minvel, str, a, hv)
  integer,                               intent(in) :: i, J
  real, dimension(NXMEM_,NYMEMQ_,NZ_),   intent(in) :: vm
  real, dimension(NXMEM_,NYMEM_,NZ_),    intent(in) :: hin
  real,                                  intent(in) :: dt
  type(ocean_grid_type),                 intent(in) :: G
  type(PointAccel_CS),                   pointer    :: CS
  real,                                  intent(in) :: maxvel, minvel
  real, optional,                        intent(in) :: str
  real, dimension(NXMEM_,NZ_), optional, intent(in) :: a, hv

! This subroutine writes to an output file all of the accelerations
! that have been applied to a column of meridional velocities over
! the previous timestep.  This subroutine is called from vertvisc.

! Arguments: i - The zonal index of the column to be documented.
!  (in)      J - The meridional index of the column to be documented.
!  (in)      vm - The new meridional velocity, in m s-1.
!  (in)      hin - The layer thickness, in m.
!  (in)      dt - The model's dynamics time step.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 PointAccel_init.
!  (in)      str - The surface wind stress integrated over a time
!                  step, in m2 s-1.
!  (in)      a - The layer coupling coefficients from vertvisc, m.
!  (in)      hv - The layer thicknesses at velocity grid points, from
!                 vertvisc, in m.

  real, pointer, dimension(:,:,:) :: &
    u_prev, &     ! The previous zonal velocity, in m s-1.
    v_prev        ! The previous meridional velocity, in m s-1.
  real    :: f_eff, CFL
  real    :: truncvel, dv
  real    :: Inorm(SZK_(G))
  real    :: e(SZK_(G)+1)
  integer :: yr, mo, day, hr, minute, sec, yearday
  integer :: k, m, mp, ks, ke
  integer :: nz
  logical :: do_k(SZK_(G)+1)
  integer :: file

!  if (.not.associated(CS)) return
  nz = G%ke
  if (CS%cols_written < CS%max_writes) then
    CS%cols_written = CS%cols_written + 1

    ks = 1 ; ke = nz ; m = 1
    do_k(:) = .false.

  ! Open up the file for output if this is the first call.
    if (CS%v_file < 0) then
      if (len_trim(CS%v_trunc_file) < 1) return
      call open_file(CS%v_file, trim(CS%v_trunc_file), action=APPEND_FILE, &
                     form=ASCII_FILE, threading=MULTIPLE, fileset=SINGLE_FILE)
      if (CS%v_file < 0) then
        call GOLD_error(NOTE, 'Unable to open file '//trim(CS%v_trunc_file)//'.')
        return
      endif
    endif
    file = CS%v_file

    if ( .not.(associated(CS%u_prev) .and. associated(CS%v_prev)) ) then
      ! Determine the index in u of um.
      do k=1,nz
        if (CS%v(i,J,k,1) /= CS%v(i,J,k,2)) then
          if (vm(i,J,k) == CS%v(i,J,k,1)) then ; m = 1
          elseif (vm(i,J,k) == CS%v(i,J,k,2)) then ; m = 2
          else
            write(file,'("Mismatched v velocities in GOLD_pointaccel")') ; exit
          endif
        endif
      enddo
      mp = 3 - m
    endif

    if (associated(CS%u_prev)) then
      u_prev => CS%u_prev
    else
      u_prev => CS%u(:,:,:,mp)
    endif
    if (associated(CS%v_prev)) then
      v_prev => CS%v_prev
    else
      v_prev => CS%v(:,:,:,mp)
    endif

    do k=1,nz
      if (((max(CS%v_av(i,J,k), vm(i,J,k)) >= maxvel) .or. &
           (min(CS%v_av(i,J,k), vm(i,J,k)) <= minvel)) .and. &
          ((hin(i,j,k) + hin(i,j+1,k)) > 3.0*G%Angstrom)) exit
    enddo
    ks = k
    do k=nz,1,-1
      if (((max(CS%v_av(i,J,k), vm(i,J,k)) >= maxvel) .or. &
           (min(CS%v_av(i,J,k), vm(i,J,k)) <= minvel)) .and. &
          ((hin(i,j,k) + hin(i,j+1,k)) > 3.0*G%Angstrom)) exit
    enddo
    ke = k
    if (ke < ks) then
      ks = 1; ke = nz; write(file,'("V: Unable to set ks & ke.")')
    endif

    call get_date(CS%Time, yr, mo, day, hr, minute, sec)
    call get_time((CS%Time - set_date(yr, 1, 1, 0, 0, 0)), sec, yearday)
    write (file,'(/,"--------------------------")')
    write (file,'(/,"Time ",i5,i4,F6.2," V-velocity violation at ",I4,": ",2(I3), &
        & " (",F7.2," E ",F7.2," N) Layers ",I3," to ",I3,". dt = ",1PG10.4)') &
        yr, yearday, (REAL(sec)/3600.0), pe_here(), i, J, &
        G%geolonv(i,J), G%geolatv(i,J), ks, ke, dt

    if (ks <= G%nk_rho_varies) ks = 1
    do k=ks,ke
      if ((hin(i,j,k) + hin(i,j+1,k)) > 3.0*G%Angstrom) do_k(k) = .true.
    enddo

    write(file,'(/,"Layers:",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(I10," ",$)') (k); enddo
    write(file,'(/,"v(m):  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (vm(i,J,k)); enddo

    write(file,'(/,"v(mp): ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (v_prev(i,J,k)); enddo

    write(file,'(/,"v(3):  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (CS%v_av(i,J,k)); enddo
    write(file,'(/,"CFL v: ",$)')
    do k=ks,ke ; if (do_k(k)) then
      CFL = abs(vm(i,J,k)) * dt * G%dx_v(i,J)
      if (vm(i,J,k) < 0.0) then ; CFL = CFL * G%IDXDYh(i,j+1)
      else ; CFL = CFL * G%IDXDYh(i,j) ; endif
      write(file,'(ES10.3," ",$)') CFL
    endif ; enddo
    write(file,'(/,"CFL0 v:",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    abs(vm(i,J,k)) * dt * G%IDYv(i,J) ; enddo

    write(file,'(/,"dv:    ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    ((vm(i,J,k)-v_prev(i,J,k))); enddo

    write(file,'(/,"CAv:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (dt*CS%CAv(i,J,k)); enddo

    write(file,'(/,"PFv:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (dt*CS%PFv(i,J,k)); enddo

    if (ASSOCIATED(CS%diag%PFv_bc)) then
      write(file,'(/,"PFvbc: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (dt*CS%diag%PFv_bc(i,J,k)); enddo
    endif
    write(file,'(/,"diffv: ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (dt*CS%diffv(i,J,k)); enddo

    if (ASSOCIATED(CS%diag%Ah_h)) then
      write(file,'(/,"Ah_h+:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%Ah_h(i,j,k)); enddo

      write(file,'(/,"Ah_h-:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%Ah_h(i,j-1,k)); enddo
    endif
    if (ASSOCIATED(CS%diag%Ah_q)) then
      write(file,'(/,"Ah_q+:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%Ah_q(I+1,J,k)); enddo

      write(file,'(/,"Ah_q-:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%Ah_q(I,J,k)); enddo
    endif

    if (ASSOCIATED(CS%diag%Kh_h)) then
      write(file,'(/,"Kh_h+:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%Kh_h(i,j,k)); enddo

      write(file,'(/,"Kh_h-:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%Kh_h(i,j-1,k)); enddo
    endif
    if (ASSOCIATED(CS%diag%Kh_q)) then
      write(file,'(/,"Kh_q+:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%Kh_q(I+1,J,k)); enddo

      write(file,'(/,"Kh_q-:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%Kh_q(I,J,k)); enddo
    endif

    if (ASSOCIATED(CS%diag%gradKEv)) then
      write(file,'(/,"KEv:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (dt*CS%diag%gradKEv(i,J,k)); enddo
    endif
    if (ASSOCIATED(CS%diag%rv_x_u)) then
      write(file,'(/,"Corv:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                 dt*(CS%CAv(i,J,k)-CS%diag%rv_x_u(i,J,k)); enddo
    endif
    if (ASSOCIATED(CS%diag%dv_dt_visc)) then
      write(file,'(/,"vbv:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
          (vm(i,J,k)-dt*CS%diag%dv_dt_visc(i,J,k)); enddo

      write(file,'(/,"dvv:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (dt*CS%diag%dv_dt_visc(i,J,k)); enddo
    endif
    if (ASSOCIATED(CS%diag%dv_adj)) then
      write(file,'(/,"dvadj: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%dv_adj(i,J,k)); enddo
    endif
    if (ASSOCIATED(CS%diag%dv_adj2)) then
      write(file,'(/,"dvadj2:",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (CS%diag%dv_adj2(i,J,k)); enddo
    endif
    if (present(a)) then
      write(file,'(/,"a:     ",$)')
      do k=ks,ke+1 ; if (do_k(k)) write(file,'(ES10.3," ",$)') a(i,k); enddo
    endif
    if (present(hv)) then
      write(file,'(/,"hvel:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') hv(i,k); enddo
    endif
    write(file,'(/,"Stress:  ",ES10.3)') str

    if (CS%split) then
      write(file,'("dvbt:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (dt*CS%v_accel_bt(i,J,k)) ; enddo
      write(file,'(/)')
      if (ASSOCIATED(CS%diag%PFv_bt)) &
        write(file,'("PFvbt: ",ES10.3)') (dt*CS%diag%PFv_bt(i,J))
      if (ASSOCIATED(CS%diag%Corv_bt)) &
        write(file,'("Corvbt:",ES10.3)') (dt*CS%diag%Corv_bt(i,J))
      if (ASSOCIATED(CS%diag%Nonlnv_bt)) &
        write(file,'("NLvbt: ",ES10.3)') (dt*CS%diag%Nonlnv_bt(i,J))
    endif

    write(file,'("h--:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') hin(i-1,j,k); enddo
    write(file,'(/,"h0-:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') hin(i,j,k); enddo
    write(file,'(/,"h+-:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') hin(i+1,j,k); enddo
    write(file,'(/,"h-+:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') hin(i-1,j+1,k); enddo
    write(file,'(/,"h0+:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') hin(i,j+1,k); enddo
    write(file,'(/,"h++:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') hin(i+1,j+1,k); enddo

    e(nz+1) = -G%D(i,j)
    do k=nz,1,-1 ; e(K) = e(K+1) + hin(i,j,k); enddo
    write(file,'(/,"e-:    ",$)')
    write(file,'(ES10.3," ",$)') e(ks)
    do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ",$)') e(K); enddo

    e(nz+1) = -G%D(i,j+1)
    do k=nz,1,-1 ; e(K) = e(K+1) + hin(i,j+1,k) ; enddo
    write(file,'(/,"e+:    ",$)')
    write(file,'(ES10.3," ",$)') e(ks)
    do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ",$)') e(K); enddo
    if (ASSOCIATED(CS%T)) then
      write(file,'(/,"T-:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') CS%T(i,j,k); enddo
      write(file,'(/,"T+:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') CS%T(i,j+1,k); enddo
    endif
    if (ASSOCIATED(CS%S)) then
      write(file,'(/,"S-:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') CS%S(i,j,k); enddo
      write(file,'(/,"S+:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') CS%S(i,j+1,k); enddo
    endif

    write(file,'(/,"u--:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') u_prev(I-1,j,k); enddo
    write(file,'(/,"u-+:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') u_prev(I-1,j+1,k); enddo
    write(file,'(/,"u+-:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') u_prev(I,j,k); enddo
    write(file,'(/,"u++:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') u_prev(I,j+1,k); enddo

    write(file,'(/,"uh--:  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (CS%uh(I-1,j,k)*G%IDYu(I-1,j)); enddo
    if (ASSOCIATED(CS%diag%uh_max)) then
      write(file,'(/," uh^--: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (CS%diag%uh_max(I-1,j,k)*G%IDYu(I-1,j)); enddo
    endif
    if (ASSOCIATED(CS%diag%uh_min)) then
      write(file,'(/," uh_--: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (CS%diag%uh_min(I-1,j,k)*G%IDYu(I-1,j)); enddo
    endif
    if (ASSOCIATED(CS%diag%uh_lay)) then
      write(file,'(/," uhL--: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (CS%diag%uh_lay(I-1,j,k)*G%IDYu(I-1,j)); enddo
    endif
    if (ASSOCIATED(CS%diag%uh_cent)) then
      write(file,'(/," uhC--: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                   (CS%diag%uh_cent(I-1,j,k)*G%IDYu(I-1,j)); enddo
    endif
    write(file,'(/," uhCp--:",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                        (0.5*u_prev(I-1,j,k)*(hin(i-1,j,k) + hin(i,j,k))); enddo

    write(file,'(/,"uh-+:  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (CS%uh(I-1,j+1,k)*G%IDYu(I-1,j+1)); enddo
    if (ASSOCIATED(CS%diag%uh_max)) then
      write(file,'(/," uh^-+: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                (CS%diag%uh_max(I-1,j+1,k)*G%IDYu(I-1,j+1)); enddo
    endif
    if (ASSOCIATED(CS%diag%uh_min)) then
      write(file,'(/," uh_-+: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                (CS%diag%uh_min(I-1,j+1,k)*G%IDYu(I-1,j+1)); enddo
    endif
    if (ASSOCIATED(CS%diag%uh_lay)) then
      write(file,'(/," uhL-+: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                (CS%diag%uh_lay(I-1,j+1,k)*G%IDYu(I-1,j+1)); enddo
    endif
    if (ASSOCIATED(CS%diag%uh_cent)) then
      write(file,'(/," uhC-+: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                              (CS%diag%uh_cent(I-1,j+1,k)*G%IDYu(I-1,j+1)); enddo
    endif
    write(file,'(/," uhCp-+:",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                  (0.5*u_prev(I-1,j+1,k)*(hin(i-1,j+1,k) + hin(i,j+1,k))); enddo

    write(file,'(/,"uh+-:  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (CS%uh(I,j,k)*G%IDYu(I,j)); enddo
    if (ASSOCIATED(CS%diag%uh_max)) then
      write(file,'(/," uh^+-: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (CS%diag%uh_max(I,j,k)*G%IDYu(I,j)); enddo
    endif
    if (ASSOCIATED(CS%diag%uh_min)) then
      write(file,'(/," uh_+-: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (CS%diag%uh_min(I,j,k)*G%IDYu(I,j)); enddo
    endif
    if (ASSOCIATED(CS%diag%uh_lay)) then
      write(file,'(/," uhL+-: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (CS%diag%uh_lay(I,j,k)*G%IDYu(I,j)); enddo
    endif
    if (ASSOCIATED(CS%diag%uh_cent)) then
      write(file,'(/," uhC+-: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (CS%diag%uh_cent(I,j,k)*G%IDYu(I,j)); enddo
    endif
    write(file,'(/," uhCp+-:",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                          (0.5*u_prev(I,j,k)*(hin(i,j,k) + hin(i+1,j,k))); enddo

    write(file,'(/,"uh++:  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (CS%uh(I,j+1,k)*G%IDYu(I,j+1)); enddo
    if (ASSOCIATED(CS%diag%uh_max)) then
      write(file,'(/," uh^++: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (CS%diag%uh_max(I,j+1,k)*G%IDYu(I,j+1)); enddo
    endif
    if (ASSOCIATED(CS%diag%uh_min)) then
      write(file,'(/," uh_++: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (CS%diag%uh_min(I,j+1,k)*G%IDYu(I,j+1)); enddo
    endif
    if (ASSOCIATED(CS%diag%uh_lay)) then
      write(file,'(/," uhL++: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (CS%diag%uh_lay(I,j+1,k)*G%IDYu(I,j+1)); enddo
    endif
    if (ASSOCIATED(CS%diag%uh_cent)) then
      write(file,'(/," uhC++: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                  (CS%diag%uh_cent(I,j+1,k)*G%IDYu(I,j+1)); enddo
    endif
    write(file,'(/," uhCp++:",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                    (0.5*u_prev(I,j+1,k)*(hin(i,j+1,k) + hin(i+1,j+1,k))); enddo

    if (ASSOCIATED(CS%diag%q)) then
      write(file,'(/,"u-eff: ",$)')
      f_eff = G%F(I-1,J)
      if (abs(f_eff) < 1e-10) f_eff = 1e-10
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (0.5*(CS%diag%q(I-1,J,k)* &
        (CS%uh(I-1,j,k) + CS%uh(I-1,j+1,k))) * G%IDYv(i,J) / f_eff); enddo

      write(file,'(/,"u+eff: ",$)')
      f_eff = G%F(I,J)
      if (abs(f_eff) < 1e-10) f_eff = 1e-10
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (0.5*(CS%diag%q(I,J,k)* &
        (CS%uh(I,j,k) + CS%uh(I,j+1,k))) * G%IDYv(i,J) / f_eff); enddo
    endif

    write(file,'(/,"D:     ",2(ES10.3))') G%D(i,j),G%D(i,j+1)

  !  From here on, the normalized accelerations are written.

    do k=ks,ke
      dv = vm(i,J,k)-v_prev(i,J,k)
      if (abs(dv) < 1.0e-6) dv = 1.0e-6
      Inorm(k) = 1.0 / dv
    enddo

    write(file,'(2/,"Norm:  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') (1.0/Inorm(k)); enddo
    write(file,'(/,"dv:    ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                    ((vm(i,J,k)-v_prev(i,J,k))*Inorm(k)); enddo
    write(file,'(/,"CAv:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                    (dt*CS%CAv(i,J,k)*Inorm(k)); enddo
    write(file,'(/,"PFv:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                    (dt*CS%PFv(i,J,k)*Inorm(k)); enddo
    if (ASSOCIATED(CS%diag%PFv_bc)) then
      write(file,'(/,"PFvbc: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                      (dt*CS%diag%PFv_bc(i,J,k)*Inorm(k)); enddo
    endif
    write(file,'(/,"diffv: ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                    (dt*CS%diffv(i,J,k)*Inorm(k)); enddo

    if (ASSOCIATED(CS%diag%gradKEu)) then
      write(file,'(/,"KEv:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                    (dt*CS%diag%gradKEv(i,J,k)*Inorm(k)); enddo
    endif
    if (ASSOCIATED(CS%diag%rv_x_u)) then
      write(file,'(/,"Corv:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
          dt*(CS%CAv(i,J,k)-CS%diag%rv_x_u(i,J,k))*Inorm(k); enddo
    endif
    if (ASSOCIATED(CS%diag%dv_dt_visc)) then
      write(file,'(/,"dvv:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
          (dt*CS%diag%dv_dt_visc(i,J,k)*Inorm(k)); enddo
    endif
    if (ASSOCIATED(CS%diag%dv_adj)) then
      write(file,'(/,"dvadj: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
          (CS%diag%dv_adj(i,J,k)*Inorm(k)); enddo
    endif
    if (ASSOCIATED(CS%diag%dv_adj2)) then
      write(file,'(/,"dvadj2:",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
          (CS%diag%dv_adj2(i,J,k)*Inorm(k)); enddo
    endif
    if (CS%split) then
      write(file,'(/,"dvbt:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                      (dt*CS%v_accel_bt(i,J,k)*Inorm(k)) ; enddo
    endif

    write(file,'(2/)')

    call flush(file)
  endif

end subroutine write_v_accel

subroutine PointAccel_init(HIS, Time, G, param_file, diag, dirs, CS)
  type(ocean_internal_state), target, intent(in) :: HIS
  type(time_type), target, intent(in) :: Time
  type(ocean_grid_type),   intent(in) :: G
  type(param_file_type),   intent(in) :: param_file
  type(diag_ptrs), target, intent(inout) :: diag
  type(directories),       intent(in) :: dirs
  type(PointAccel_CS),     pointer    :: CS
! Arguments: HIS - For "GOLD Internal State" a set of pointers to the fields and
!                  accelerations that make up the ocean's physical state.
!  (in)      Time - The current model time.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      diag - A structure containing pointers to common diagnostic fields.
!  (in)      dirs - A structure containing several relevant directory paths.
!  (in/out)  CS - A pointer that is set to point to the control structure
!                 for this module
  character(len=128) :: version = '$Id$'
  character(len=128) :: tagname = '$Name$'
  character(len=40)  :: mod = "GOLD_PointAccel" ! This module's name.

  if (associated(CS)) return
  allocate(CS)

  CS%diag => diag ; CS%Time => Time

  CS%u => HIS%u ; CS%v => HIS%v ; CS%uh => HIS%uh ; CS%vh => HIS%vh
  CS%CAu => HIS%CAu ; CS%CAv => HIS%CAv ; CS%PFu => HIS%PFu ; CS%PFv => HIS%PFv
  CS%diffu => HIS%diffu ; CS%diffv => HIS%diffv
  CS%T => HIS%T ; CS%S => HIS%S ; CS%pbce => HIS%pbce
  CS%u_accel_bt => HIS%u_accel_bt ; CS%v_accel_bt => HIS%v_accel_bt
  CS%u_prev => HIS%u_prev ; CS%v_prev => HIS%v_prev
  CS%u_av => HIS%u_av; if (.not.associated(HIS%u_av)) CS%u_av => HIS%u(:,:,:,1)
  CS%v_av => HIS%v_av; if (.not.associated(HIS%v_av)) CS%v_av => HIS%v(:,:,:,1)

  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mod, version, tagname, "")
  call get_param(param_file, mod, "U_TRUNC_FILE", CS%u_trunc_file, &
                 "The absolute path to the file where the accelerations \n"//&
                 "leading to zonal velocity truncations are written. \n"//&
                 "Leave this empty for efficiency if this diagnostic is \n"//&
                 "not needed.", default="")
  call get_param(param_file, mod, "V_TRUNC_FILE", CS%v_trunc_file, &
                 "The absolute path to the file where the accelerations \n"//&
                 "leading to meridional velocity truncations are written. \n"//&
                 "Leave this empty for efficiency if this diagnostic is \n"//&
                 "not needed.", default="")
  call get_param(param_file, mod, "MAX_TRUNC_FILE_SIZE_PER_PE", CS%max_writes, &
                 "The maximum number of colums of truncations that any PE \n"//&
                 "will write out during a run.", default=50)
  call get_param(param_file, mod, "SPLIT", CS%split, &
                 "Use the split time stepping if true.", default=.true.)

  if (len_trim(dirs%output_directory) > 0) then
    if (len_trim(CS%u_trunc_file) > 0) &
      CS%u_trunc_file = trim(dirs%output_directory)//trim(CS%u_trunc_file)
    if (len_trim(CS%v_trunc_file) > 0) &
      CS%v_trunc_file = trim(dirs%output_directory)//trim(CS%v_trunc_file)
    call log_param(param_file, mod, "output_dir/U_TRUNC_FILE", CS%u_trunc_file)
    call log_param(param_file, mod, "output_dir/V_TRUNC_FILE", CS%v_trunc_file)
  endif
  CS%u_file = -1 ; CS%v_file = -1

end subroutine PointAccel_init
end module GOLD_PointAccel
